Core reconstruction in pseudopotential calculations 
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A new method is presented for obtaining all-electron results from a pseudopotential calculation. 
This is achieved by carrying out a localised calculation in the region of an atomic nucleus using the 
embedding potential method of Inglesfield [J.Phys. C 14, 3795 (1981)]. In this method the core 
region is reconstructed, and none of the simplifying approximations (such as spherical symmetry 
of the charge density/potential or frozen core electrons) that previous solutions to this problem 
have required are made. The embedding method requires an accurate real space Green function, 
and an analysis of the errors introduced in constructing this from a set of numerical eigenstates 
is given. Results are presented for an all-electron reconstruction of bulk aluminium, for both the 
charge density and the density of states. 

PACS numbers: 71.15.Hx, 71.15.Ap, 71.15.-m 



I. INTRODUCTION 

Total energy pseudopotential methods have taken 
pride of place in the first principles simulation of con- 
densed matter in recent years due to their efficient use of 
computing resources and their suitability for structural 
optimisation However, the charge density resulting 
from a pseudopotential calculation is incorrect in the re- 
gion near the atomic nuclei - it does not include core 
states, and the valence states have the wrong structure. 
In order to obtain accurate values for any quantity that 
depends on the true charge density, such as hyperfme in- 
teractions or X-ray structure factors, we must obtain the 
correct electron charge density from the pseudopotential 
calculation. Other methods are available that calculate 
the states of all the electrons in the system (Full-potential 
Linear Augmented Plane Wave (FLAPW)[2j, Projector 
Augmented Plane Wave (PAW)|H, Linear Muffin Tin 
Orbital (LMTO) 3 , KKR Green function [| and Tight 
Binding Q), but these methods tend to be more com- 
putationally expensive, not as well suited to structural 
optimisation, or are less accurate than pseudopotential 
methods. 

In view of this it is desirable to extend the pseudopo- 
tential method by adding an extra step after the pseudo- 
system has been solved, ie to choose an atom for which 
we require the core and correct valence states, and re- 
construct these correct states. It would be hoped that 
solving for one atom with different boundary conditions 
on a sphere surrounding it would be fairly straightfor- 
ward. However, a number of difficulties quickly present 
themselves. The purpose of this paper is to describe 
and validate a new procedure for carrying out this core 
reconstruction that makes essentially the same physical 
approximations as the FLAPW method, and so can be 
expected to provide the same accuracy. 

This problem has been addressed by several previous 
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workers. In their paper Gardner and Holzworth [tJ re- 
construct the correct states for isolated Si and Ru atoms 
from the pseudo-atom, by applying direct integration, 
effectively 'inverting' the pseudopotential. They obtain 
good results, showing that this reconstruction approach 
is at least possible. However, reconstructing the states of 
an atom in a lattice is considerably more complex. For 
this case the valence states form a continuum, so the re- 
construction must be fitted into the band structure of 
the lattice, and in addition to this the potential is not 
spherically symmetric. Vackaf and Simunek Q describe 
a method for reconstructing the states for a pseudo-atom 
within a lattice. Their method relies on direct integra- 
tion and assumes the charge density, boundary conditions 
and self-consistent potential are spherically symmetric, 
although the core states are allowed to relax. The er- 
rors in the resulting eigenfunctions are fairly large, al- 
though they do obtain the correct nodal form for the 
eigenstates. Kuzmiak et alQ perform a pseudopoten- 
tial calculation, and orthogonalise the resultant pseudo- 
states to the original core states. This would work for the 
original formulation of pseudopotential methods, where 
the pseudopotential is defined in this way, but for modern 
norm-conserving pseudopotentials this does not give the 
correct solution to the problem and the errors present are 
difficult to control or even quantify. The most complete 
solution to the problem presented so far is that due to 
Meyer et al [lfj. In their method the correct states are 
reconstructed by direct integration. In order to decouple 
the radial wave equations for the reconstruction calcu- 
lation they make the assumption of spherical symmetry 
of the self consistent potential, but asymmetric bound- 
ary conditions for the valence states are allowed. Within 
their scheme the core is still frozen. 

In this paper a new method for performing this kind of 
core reconstruction is described that does not make any 
of the assumptions of previous approaches. The method 
presented here follows a different path to achieving the re- 
construction, does not require spherical averaging of the 
self-consistent potential, provides an aspherical charge 
density, and does not assume a frozen core. The first 
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step is to carry out a plane-wave pseudopotential calcu- 
lation for the system of interest, and construct the single 
particle Green function for the valence electrons present 
in this system. This Green function is then used to create 
an embedding potential, as described by Inglcsficld 11]. 
An all-electron localised atomic calculation is then car- 
ried out in a space containing the core region of the atom 
of interest, with the embedding potential taking into ac- 
count the rest of the atoms in the lattice. Effectively the 
system that is solved for is one all-electron atom in a lat- 
tice of pseudo-atoms. This approach preserves the full 
generality and flexibility of the pseudopotential method. 
In addition, the all-electron calculation is carried out only 
for the atom(s) of interest, so computational effort can 
be applied sparingly. 

In the next section a brief description is given of the 
embedding approach, how it is applied in this case, and 
how it relies on an accurate knowledge of the Green func- 
tion for the substrate system. Appendix [A] gives a more 
complete description of the method and its properties. 
In section Hm we describe how the spectral representation 
can be efficiently applied in order to construct an accu- 
rate Green function, and how the incompleteness of the 
available set of states introduces a significant error that 
must be corrected. Combining the results of these two 
sections we then obtain the embedding potential from the 
Green function and perform a localised calculation using 
the embedding terms to take into account the rest of the 
lattice, as described in section IIVI Section [V] gives the 
results of a reconstruction carried out for bulk FCC alu- 
minium, and the convergence of the method is discussed. 
Rydberg atomic units are used throughout the paper. 



II. EMBEDDING 

An embedding method can be thought of as solving 
for a system in a sub-domain of space, denoted region / 
in what follows, where the influence of the system out- 
side of this sub-domain (region II) is taken into account 
from a previous solution, and is not recalculated. This is 
shown in Fig. [1] for the core reconstruction case consid- 
ered here. The embedding surface, S, separates regions 
I and 77. Region I is a sphere of radius r s enclosing 
the core region of the atom where the pseudo-states are 
incorrect and region // is the remainder of the lattice of 
pseudo-atoms. It is implicitly assumed throughout this 
work that norm-conserving pseudopotentials are used, so 
that wavefunctions between core regions are correct. The 
embedding surface S is assumed to be in such a region; 
it follows that core regions do not overlap, hence r s > r c 
where r c is the largest pseudopotential core radius. 

Using Inglesfield's method ll| an 'embedding poten- 
tial' is obtained from the substrate system, and added to 
the Hamiltonian for the embedded region. This embed- 
ding potential ensures that the states of the system in 
the embedded region satisfy the correct boundary con- 
ditions. Inglesfield's method has the advantage that it 




FIG. 1: Geometry of embedding calculation. Region II is the 
substrate region, region I is the embedded region, and the 
all-electron states are reconstructed from a knowledge of the 
pseudo-states on the surface 5. 



requires knowledge of the properties of the substrate sys- 
tem only on the surface separating the embedded and 
substrate regions. In appendix [A] a brief description of 
the derivation of the embedding method is given, followed 
by two expressions for the embedding potential in terms 
of the Green function for the substrate system, one of 
which has not appeared in the literature before. 

Within a continuum the Hamiltonian of the embedded 
system takes the form 



H emb (E) =Hj + 6(r s - r) 



^--T(r s ,r' s ;E) 
an s 



(1) 



where H em b(E) is the embedded Hamiltonian that yields 
the states with correct boundary conditions, Hi is the 
normal Hamiltonian for the embedded region, and T is 
the embedding potential. It should be understood at this 
point that T(r s ,r' s ;E) acting on a function denotes the 
integration over the surface S, as described in appendix 
IA1 Equation ([1]), when solved in region I, will give the 
correct solution for the system represented by regions 
/ and II, with region II represented entirely through 
the embedding potential term. The embedding potential 
required in equation (Eq. |T])) is given by the operator 



on' 



(2) 



where Q is the matrix representation of the Green func- 
tion on the surface S in terms of a set of basis func- 
tions orthonormal over the surface, and the derivative is 
the normal derivative of Q outward from the surface and 
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with respect to the second spatial variable of the Green 
function. A derivation of Eqs. |T]) and $2$ is given in 
appendix [A] 



III. GREEN FUNCTION AND EMBEDDING 
POTENTIAL 

The embedding method of the previous section is to be 
applied with the substrate system represented by the re- 
sults of a plane-wave pseudopotential calculation, hence 
the Green function on the embedding surface must be 
obtained from the Bloch states expanded in plane-waves. 
The natural way of constructing this Green function is via 
the spectral representation, and a finding good approx- 
imation to this spectral representation is the concern of 
this section. 

For the periodic system the states are characterised by 
two quantum numbers, the discrete band index n and the 
continuous crystal momentum k, so the spectral repre- 
sentation takes the form 



„ Jbz 



E n (k)-E 



(3) 



where ^(r) is a Bloch state, E n (k) is its eigenenergy and 
E is complex in general. For the periodic lattice both the 
total number of states and the number of states within 
a given energy interval (in a band) is infinite, but any 
numerical calculation can only provide states at a finite 
number of k-points, and for a finite number of bands 
In view of these restrictions the approximation of 
the spectral representation falls naturally into two parts 
- approximating the Brillouin zone integral from a finite 
number of k-points, and approximating the infinite band 
sum. It should be noted that although it is well estab- 
lished that Green function methods and the spectral rep- 
resentation can be defined within a limited basis set [l2| , 
the embedding method cannot be applied in this way 
since it relies on the properties of the Green function in 
real space [l3j . 



A. The Spectral function 

The simplest way to apply the defining equations for 
the embedding potential is to expand the Green function 
in a set of basis functions that are orthogonal over the 
embedding surface. For the core reconstruction consid- 
ered here the embedding surface is a sphere centred on 
the atomic site of interest, hence a natural set of basis 
states are the spherical harmonics. First the Kohn-Sham 
wavefunctions are expanded on the surface S, 



g 



L 



c4 n) (k)y L (r) (4) 



where L = (Im) , the combined index of the spherical har- 
monic Yl, and C™ (k) are the expansion coefficients of the 



eigenstates in the plane-wave representation. The expan- 
sion coefficients (k) can be found using the identity 

M 



e** = 4nJ2i l MQr)Y£(q)Y L (r) 



(5) 



hence 



a[ n \k) = 47ri , $>(|k + g|r i )y£(k + g)C£(k). (6) 
g 

From these a 'spectral function' J-ll'{E) = 
ilm [Gll'(E)] can be defined, where Gll'(E) are 
the expansion coefficients for the Green function. This 
spectral function is given by the equation 

T LU {E) = V / d 3 k 4 n) (k)4^*(k)5(£ - E n (k)) 
n Jbz 

(7) 

and is related to the Green function by the convolution 
integral 



Gll'(E) 



— OG 



E' — E 



(8) 



In Eq. ([7]) the integral on the right hand side reduces to 
the surface integral 



T LL ,{E)=Y. ! 



a L ( k ) a L> ( k ) 



E=E n (k) 



|V k £| 



dS (9) 



due to the delta function, hence evaluation of the spec- 
tral representation reduces to the evaluation of a surface 
integral in k-space and a singular convolution integral in 
energy space. The surface integral is carried out usin g th e 
linear analytic tetrahedron method [HI, EE E3, ll8Ul9t , 
which results in a spectral function that has the correct 
analytic structure in that it is continuous in energy. Only 
the irreducible wedge need be sampled with the symme- 
try of the crystal used to complete the rest of the integral 

In order to evaluate Eq. ([5]) the singular integral itself 
must be approximated from a sampling of the function 
at a finite number of energy values. Here the spectral 
function is interpolated between consecutive energy val- 
ues, and the contribution to the integral from each of 
these ranges calculated. We use a polynomial interpola- 
tion, so the integrals can be evaluated analytically [2lj . 
This results in an approximation to the Green function 
that can be evaluated for complex E and which has the 
analytic properties appropriate for a continuum of states; 
it is essentially a generalisation of the approach used by 
Kuzmiak et al [9j to complex energies. An alternative 
application of the linear analytic tetrahedron method to 
the calculation of Green functions is given Lambin and 
Vigneron [12] where Eq. §5§ is evaluated for each tetrahe- 
dron analytically, within the linear interpolation scheme. 
Although this approach is more direct and introduces 
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none of the errors due to approximating Eq. ([8]), it re- 
quires the interpolation to be applied to all tetrahedra, 
whereas the surface integral in Eq. ([9|) requires only a 
sub-set of tetrahedra to be considered for each energy. 
In addition to this the surface integral and convolution 
route allows a greater flexibility in the degree of approx- 
imation applied, as discussed below. 



B. Completing the Incomplete Set of States 

So far the evaluation of the Green function takes into 
account the continuum nature of the band states, but 
for the spectral representation to describe a real space 
Green function the set of states used in Eq. ([3]) must 
be complete. We have found that in order to obtain ac- 
curate Green functions using the spectral representation, 
we must include the complete set of n s states associated 
with the plane- wave basis set (for further discussion see 
Trail [13] )• n % is the number of plane- waves in the ba- 
sis set, determined by the plane- wave energy cut-off, and 
these states will be obtained by direct matrix diagonali- 
sation of the Kohn-Sham Hamiltonian. This set of states 
is still not complete in real space, hence we must con- 
sider the errors introduced by not including the energy 
bands that are not available in a finite plane-wave cal- 
culation. In what follows we refer to the spectral repre- 
sentation that includes only the n g lowest energy bands 
as the incomplete spectral representation. It is impor- 
tant to realise that we mean the states are incomplete in 
real space, they are, of course, complete in the sub-space 
spanned by the plane-wave basis set. 

It is not immediately apparent that this incomplete- 
ness will have a significant effect, and it could be hoped 
that the 'missing' high energy states are so far above the 
energies of interest (ie at or below the Fermi energy) that 
any error introduced by their absence will be negligible. 
This is only partly true, and the properties of the error 
introduced by incompleteness are derived in appendix [B] 
for a free electron gas. Correcting for this incompleteness 
not only speeds convergence with respect to n g , it also 
ensures that the approximation of the Green function has 
the correct analytic form. We follow James and Wood- 
ley [23l ] and approximate the high energy states 'missing' 
from the incomplete spectral representation by free elec- 
tron states. The plane-wave basis set used to represent 
the lowest n s bands is described by |g| 2 < E max , where 
g are reciprocal lattice vectors and E max is the standard 
plane-wave cut-off energy. Consequently, the free elec- 
tron states required to 'top up' the incomplete spectral 
representation are those described in the reduced zone 
scheme by |g| 2 > E max and k in the first Brillouin zone. 

In order to calculate the required correction we calcu- 
late an incomplete spectral representation for free elec- 
trons with exactly the same basis as for the pseudopo- 
tential states, and subtract this from the analytic free 



electron Green function. This yields the approximation 

G « Q p E s Zt° - GeZ. + g ™ ee ( 10 ) 

where the first term on the right hand side is the in- 
complete spectral representation of the pseudo-states, the 
second is the incomplete spectral representation of free 
electron states and the final term the complete spectral 
representation for free electron states (ie the analytic free 
electron Green function) . In terms of the spectral func- 
tion, and the convolution integral used to transform 
this into the spectral representation, Eq. © becomes 

/^-pseudo 1 _ ■jrfree / p/\ 

W V-E dE> + GTl ' {E) 

(11) 

where Q is now a Green function with the correct analytic 
form, E is complex and E' is real. J^P seudo and J^J ree are 
the spectral functions associated with the pseudo-states 
and the free space states respectively, calculated with n s 
basis functions. The last term, Q°°, is the analytic free 
space Green function given by [3, H3 

&lv = ikji(kr s )hi(kr s )S w (12) 

where ji is the spherical Bessel function of the first kind, 
hi is the spherical Hankel function of the first kind and 
E = k 2 + Vq, where Vq is the average potential within the 
unit cell. In practice the specific value of Vq is not criti- 
cal since the error term varies only slowly with the Green 
function energy, as discussed in appendix [B] The normal 
derivative of the Green function can be obtained by tak- 
ing the radial derivative of the Green function, and Eq. 
@ directly applied to give the embedding potential (a 
factor of -5 must be included to normalise the spherical 

harmonics over the embedding sphere). Eq. (|A12j) wa s 
not used due to problems with small errors in the Green 
function introducing anomalous poles at the bottom of 
the lowest band. 



IV. EMBEDDING CALCULATION 

In this section we describe the all-electron calculation 
carried out in the region near the nucleus, using the 
embedding potential described above. The normal den- 
sity functional framework is used, with the Kohn-Sham 
Hamiltonian extended by the addition of the embedding 
terms. Since these extra terms are functions of energy 
the eigenvalue solution of the Hamiltonian is not simple, 
hence the charge density is obtained directly from the 
Hamiltonian via the Green function method as described 
by Williams et al [12]. The method employed in the 
all-electron calculation is similar to that used by other 
workers (eg Trioni et al 25}), but generalised so as not to 
require any particular symmetry of the charge density or 
self-consistent potential, and to include core electrons. 
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The Embedded Hamiltonian 



region A results in 



Each basis function is a product of a radial function 
and a spherical harmonic, with the radial part defined 
as an augmented plane-wave for r < s (region B) and 
a spherical Bessel function for s < r < r s (region A), 
where s is a parameter of the calculation. These are cho- 
sen since Linearised Augmented Plane Waves (LAPW) 
orbitals describe the all-electron valence states well near 
the nucleus, and in the region nearer the embedding ra- 
dius the spherical Bcsscl functions provide the flexibility 
required to satisfy the boundary conditions. 

The basis functions in region B are found by solv- 
ing the Dirac equation using numerical integration in 
the spherical part of the Kohn-Sham potential. The 
method used is that described by Koelling and Harmon 
[26| . where the Dirac equation is approximated in the 
form of a Schrodinger equation which does not include 
spin-orbit interaction but takes other relativistic effects 
into account. A solution is found at a fixed 'pivot' en- 
ergy, E p , and the resultant radial functions are denoted 
ui{r). The energy derivative of these functions are also 
obtained and orthogonalised to the associated u;(r), and 
these orthogonalised energy derivatives are denoted ui (r) 
(see Krasovskii [2Z|, or Takeda and Kiibler [H). In A 
the radial basis functions used are spherical Bessel func- 
tions of the first kind, defined as ji(gir) where gt = 22 
The index i is chosen to take integer values, to give a 
set of functions, and d is a parameter of the calculation 
(larger than r s to allow sufficient flexibility for the basis 
function to satisfy arbitrary boundary conditions). The 
radial parts of the basis functions, \u ( r ) are therefore 



a-uui(r) + b a ui(r) < r < s 
ji {gir) s <r <r s 



(13) 



Parameters an and bu are chosen to ensure xu is contin- 
uous in amplitude and derivative at the AB boundary. 
Typically r s is chosen to be half the distance between 
nearest neighbour atoms in the lattice, s ~ 0.9r s and 
d ~ 2r s to give a good description of the valence elec- 
trons. For i the range 1 ... 4 is typical. 

The embedded Hamiltonian matrix is then expanded in 
terms of these basis functions, for a Kohn-Sham potential 
given by 



V(v)=Y,V L {r)Y L {v) 



(14) 



where L = (Im) , the index of the spherical harmonic. We 
write Hamiltonian matrix, H em b, as the sum of 3 parts, 



TI II 1 



(15) 



with H A the contribution from region A, H B the con- 
tribution from region B, and £ the contribution from 
the embedding terms at the surface iS. Integrating over 



H iL,jL> 



g) / r 2 dr\ji(gir)ji(gjr)]S 



LL> 



+ V S L L , L „ I r 2 dr \ji{ 9i r)V LII (r)j v (g 3 r)} , (16) 



where S^, L „ = J Y£YL'YL"dV, are Gaunt coefficients. 
For region B we find a sum of spherical and aspherical 
terms 



H 



iLJL' 



[auct,jiEp{ui | m) 

+aubji(ui | iii) + b a bjiE p (ui \ it;)] Sll' 



r 2 dr X a(r)V L „{r)x 3 v(m 



where E p is the 'pivot' energy at which ui(r) is calculated, 
and the bra-ket's denote integration over region B only. 
The integrals in the above expressions are carried out 
analytically or numerically as appropriate. Taking the 
normal derivative and embedding potential terms in Eq. 
(fTj) we obtain the contribution from the embedding terms, 
S, as 



HL,jL' = ji{gi'r s )T LL ,{E)j v {g J r s ) 
+3j3li9i r s)j'i'(9jrs)h,L' 



(18) 



In addition to the Hamiltonian matrix the overlap ma- 
trix of the basis functions, O, is also required. This is 
given by 



(19) 



OtLjL' =5ll> I ji(g t r)ji(gjr)r dr 

J s 

+ [aua,ji(ui | ui) + bubji(iii | iii)} 8 LL > 



where the first integral on the right hand side can be 
performed analytically [H| and the second numerically. 



B. The Embedded Green function and Charge 
Density 

The Green function of the embedded system is ob- 
tained directly by matrix inversion, 



G(E) = (H emb (E) - EOy\ 



(20) 



where H„ 



is the embedded Hamiltonian described 



above. To obtain the local density of states, n(r;E), 
we apply the identity 29] 



n(r;E) = -Im G(r,r;E). 

TT 



(21) 



By filling all states to the Fermi energy (which is ob- 
tained from the original pseudopotential calculation us- 
ing the linear analytic tetrahedron method) we can ob- 
tain the valence charge density. Unfortunately the num- 
ber of points required to evaluate the charge density by 
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integrating the local density of states is fairly high, due 
to its fine structure. This can be avoided be taking ad- 
vantage of Cauchy's theorem [H, EH and the analytic 
properties of the Green function. We choose a contour C 
from some energy below the lowest eigenvalue, following 
a half circle into the positive complex plane and termi- 
nating on the real axis at the Fermi energy, so the charge 
density expanded in spherical harmonics is given by 



ijLL' 

x J c (G iLjL ,S£' L „ - G* LjL ,Si, v ) (22) 

where the GiLju are evaluated at the complex energies 
on the contour. In this equation S^, L „ is a Gaunt coef- 
ficient, and is zero for I" > 21 or 21', so a l ma x x lmax 
Hamiltonian matrix will result in a charge density con- 
taining components I < 2l max . It is due to the basis 
functions being complex that this expression does not 
simply involve the imaginary part of the Green function 
matrix. As the imaginary part of the energy increases the 
Green function becomes smoother and more featureless, 
hence a more sparse sampling is sufficient to approximate 
this integral accurately. Gaussian integration [301 ] is used 
to perform the integral along contour C, and converged 
integrals are typically obtained for around 16 points on 
the contour. 

The core states could be obtained by finding the dis- 
crete states of the embedded Hamiltonian. but this would 
be cumbersome (see appendix |A")). Core states are there- 
fore obtained using the Kohn-Sham potential within the 
embedding sphere, and a constant potential outside of the 
sphere. This gives accurate results as the core states are 
strongly localised within the core region, so the poten- 
tial outside of the embedding sphere is largely irrelevant. 
Only the spherically symmetric part of the potential is 
used to obtain the states, and the constant potential out- 
side the sphere is taken to be continuous with the spher- 
ical part. This spherically symmetric approximation has 
been used by many workers (eg Blaha et al [3l[ ; Methfes- 
sel and Frota-Pesso |32j), and can easily be extended to 
included aspherical effects perturbatively as described by 
Ehmann and Fahnle [33| . Sternheimer [34j and Lauer et 
al [35| . To obtain the core charge density a fully relativis- 
tic treatment is used, solving the Dirac equation within 
the embedding region. 

This approach allows the construction of core-states 
that are self-consistent within the embedded system, 
however the influence of core relaxation on the substrate 
system (and so on the embedding potential) is not consid- 
ered. Essentially we are assuming that the original pseu- 
dopotential approximation is valid, and that the valence 
pseudo-states are accurate outside of the core region. In 
order to go beyond this assumption it should be possi- 
ble to obtain a new pseudopotential from the embedded 
core states and to iterate to full self-consistency, but this 
has not been implemented. Similar considerations are 



discussed by Blochl It should be noted that for the 
embedding method applied here the core states are truly 
localised, since we are embedding an all-electron atom at 
one site within a lattice of pseudo-atoms. 



C. Self Consistency 

The Kohn-Sham potential, V K (r), is calculated from 
the total charge density via the Local Density Approxi- 
mation (LDA), and takes the form 



(r) = V? uc {r) + V? art (r) + V* c {r), (23) 



a sum of the nuclear potential, the Hartree potential, 
and the exchange-correlation potential respectively. The 
nuclear potential is assumed to be that of a point charge 
centred at the nucleus, and the Hartree potential, 



Hart 



47T 

+27TT' 

+a L r l , 



r' 2 dr' 
r dr 



(24) 



is obtained from Poisson's equation. The constants at 
are defined by the boundary condition that the Kohn- 
Sham potential, Vj^ s (r), be equal to the self-consistent 
potential of the original pseudopotential calculation on 
the embedding sphere. 

The LDA is used to obtain the exchange-correlation 
potential from the charge density using the Perdew and 
Zunger parameterisation (3(| of the results of Ceperly 
and Alder [37| (this was also used in the original pseu- 
dopotential calculation). V xc (r) is a non-linear func- 
tion of the charge density, so V* (r) cannot be directly 
obtained from Pi(r). However, deviations from a spher- 
ical density are expected to be small, so we expand the 
exchange-correlation potential as a Taylor series 



V xc (p) = V^( P0 ) + 



XC I 



gyXC 



dp 



(p - Pa) 



(25) 



with po the spherical part of the charge density. A 
truncation to linear order is sufficiently accurate (the 
quadratic terms are negligible for the systems considered 
to date) and can be applied directly. Convergence to self- 
consistency is achieved, with instabilities controlled and 
convergence speed increased by applying Broyden mixing 



V. RESULTS 

In order to test the effectiveness of the embedding 
method for core reconstruction we examine aluminium 
in the FCC structure. We begin by giving an example 
of the calculated embedding potentials, and some con- 
crete evidence that the approximations made can give an 
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accurate embedding potential. The initial self-consistent 
plane-wave calculation is carried using the LDA, a plane- 
wave cut-off of 400 eV and 60 k-points within the FCC 
irreducible wedge. These parameters are more than ad- 
equate to obtain effectively perfect convergence, so we 
can attribute any errors in the reconstruction to the 
embedding method. A Kerker [40j pseudopotential is 
used with a core radius r c — 2.19 au, small enough 
for the core regions not to overlap. The resulting self- 
consistent potential is then used to obtain the full set of 
eigenstates by matrix diagonalisation for a plane-wave 
basis set |g| 2 < E max , and at the k-points required 
to carry out the Brillouin zone integral using the lin- 
ear analytic tetrahedron method. These states are then 
used to construct the embedding potential, with the em- 
bedded sphere taken as the 'touching sphere' radius, 
r s = 2.705 au. The convergence behaviour of the em- 
bedding potential depends on the number of plane- waves 
in the basis, n g , the number of k-points in the irreducible 
wedge of the FCC Brillouin zone, n^, and the spacing of 
energy points (AE) used to sample the spectral function, 
T (see Eq. ©). 

Fig. [5] shows the r( 10 10 ) matrix element for the imagi- 
nary part of the energy equal to 0.1 eV. Results are shown 
for E max — 200 eV and AE — 0.3 eV with linear interpo- 
lation of the spectral function, and for both 240 and 505 
k-points in the irreducible wedge of the Brillouin zone. 
There is very little discernable difference between the two 
results, < 0.002 for both the real and imaginary parts, 
hence we consider 505 k-points as effectively converged. 
Using 89 k-points results in larger errors (< 0.01), but 
was found to be adequate for performing an accurate 
reconstruction, as discussed below. The difference be- 
tween E max — 200 and 400 eV is even smaller (< 0.001), 
and comparing AE — 0.3 eV to 0.1 eV again gives a 
difference of < 0.001. The other matrix elements show a 
similar convergence behaviour. Embedding potentials for 
use in the reconstructions were therefore calculated with 
E max = 200 eV, AE = 0.3 eV and 240 k-points, for en- 
ergies required on the contour of integration (Eq. (|22l) ~). 
Reconstructions were performed for l max = 6, d = 4.0, 
4 Bessel functions in the basis (Eq. (TT3|) ), and 16 points 
along the contour of integration. These provided the to- 
tal charge density within the radius r s , the density of 
states within the embedded region, the self consistent 
potential and the core eigenstates. 



A. Embedding with a Pseudopotential 

Before reconstructing the all-electron charge density, 
we carry out a reconstruction in which no core states 
are included and the core and nucleus are described by a 
pseudopotential, the same Kerker pseudopotential used 
in the original plane-wave calculations. This provides a 
stringent test of the accuracy and reliability of the entire 
embedding approach, since for a successful implementa- 
tion the resulting valence charge density should be the 





range peak 


(xl(T 4 e an' 3 ) 


R(%) 


Pseudo 


r < r s 


-5.56 


0.49 


Pseudo 


r c < r < r s 


-5.56 


0.48 


All-electron 


r c < r < r s 


-5.47 


0.46 



TABLE I: Errors in charge density for each reconstruction. 
Pseudo refers to the reconstruction performed using a pseu- 
dopotential to represent core states, and All-electron refers to 
the full all-electron reconstruction. The pseudopotential core 
radius is r c — 2.19 au, and the embedding sphere radius is 
r a — 2.705 au. 



same as the original plane- wave charge density through- 
out the embedding sphere. The Fermi energy used within 
the reconstruction is obtained from the original plane- 
wave calculation, hence the total valence charge within 
the embedding sphere will be different from the plane- 
wave calculation, and the size of this difference provides 
a first indication of how accurate the reconstruction is. 
The original plane- wave calculation has 2.298 electrons 
within the embedding sphere, whereas the reconstruc- 
tion gives 2.300 - an error of 0.08 %. This close agree- 
ment suggests a successful reconstruction, however this 
is a fairly gross measure of success since it takes no ac- 
count of the structure of the electron states within the 
embedding sphere and compares only the spherical part 
of the charge density. 

A more demanding measure of the accuracy of the re- 
construction is to compare the charge density and/or self 
consistent potential of the original pseudopotential cal- 
culation with the reconstruction. The error in the charge 
density is quantified as the peak error and the R factor 
[30l | defined as 

R = J \p recon (r)-pP seud {r)\d 3 r/ J \p pseud (v)\a* v , (26) 

and the errors over the whole embedding region are given 
in Tablefl] We also give the error within the shell between 
the core radius of the pseudopotential and the embedding 
sphere, which will allow comparison with the all-electron 
results. 

Fig. shows a contour plot of p(r), for both the orig- 
inal plane-wave and reconstructed charge density. These 
are taken in the {100} plane, with the embedded atom at 
the centre and the square shown just enclosing the em- 
bedding sphere. There is excellent agreement between 
the original and reconstructed densities, again demon- 
strating a successful reconstruction. 



B. The Reconstructed All-electron Charge Density 

The total valence charge within the embedding sphere 
is 2.299 electrons for the all-electron reconstruction, 
which again agrees very well with the original plane- wave 
result. As above we also compare the charge density of 
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FIG. 2: (a) Real part, (b) imaginary part of aluminium embedding potential matrix element, T^o.io), for — 240 (solid line, 
for clarity shown offset by —0.02 for both real and imaginary parts) and nk = 505 (dashed line). The difference is also shown 
as the solid line near rV 10 ,io) = 0-0- 



the original pseudopotential calculation with the recon- 
struction. For a successful reconstruction, the original 
pseudopotential charge density should agree with the re- 
constructed charge density for r c < r < r s . Fig. \3jp 
shows p val (r) along a line in the [011] direction for both 
the original plane-wave and all-electron reconstruction 
calculations, with the reconstructed atom at the origin. 
Agreement between the two results is excellent in the re- 
gion r > r c . The errors between r c and r s are given in 
Table QJ and these are similar to those of the pseudopo- 
tential reconstruction over the same region. 

As well as the charge density, errors in the self- 
consistent potential have also been considered. The 
reconstructed potential should agree with the original 
(plane-wave) potential outside of the core radius. The R- 
factor is smaller than for the charge density (~ 0.06%), 
largely due to the fact that the reconstructed and orig- 
inal self consistent potential are forced to agree at the 
embedding sphere. 

C. Original Pseudo and Reconstructed Density of 
States 

Another quantity whose accuracy is important is the 
Density of States (DOS). Since there is a 1 : 1 correspon- 
dence between the all-electron valence and pseudo-states, 
the eigenvalues are equal and the pseudopotential is norm 
conserving, then a successful reconstruction should yield 
an identical DOS to the original results. 



The DOS within r s is calculated from the plane- wave 
results by applying the tetrahedron method to obtain the 
local density of states, and integrating this over the vol- 
ume contained by the embedding sphere. For the recon- 
structed DOS a self consistent calculation is carried out, 
and the embedded Green function obtained along a con- 
tour parallel to the real axis, with an imaginary energy 
of 0.1 eV. The LDOS is taken from the imaginary part 
of this Green function, and integrated over the volume of 
the embedding sphere, hence the reconstructed DOS is 
smoothed by a Lorentzian of width 0.1 eV, small enough 
for the fine structure to be apparent. In Fig. [4] the DOS 
within the embedding sphere is shown, for pseudo and 
all-electron reconstructed states. Agreement is excellent, 
with a maximum error of ~ 1 %. 



D. Convergence 

The results given above for the all-electron reconstruc- 
tion were reproduced for a number of different parame- 
ters of the embedding potential and reconstruction cal- 
culation to assess the convergence of the reconstructed 
results. In what follows we define a charge density as 
'well converged' when the associated R value is less than 
0.5 %. 

First we consider the reconstruction calculation itself. 
Convergence was easily achieved with respect to the pa- 
rameters of the reconstruction calculation, that is, the 
number of basis functions in Eq. (|13p and the number 
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FIG. 3: (a) Contours of constant valence charge density in the {100} plane for the original plane- wave (dashed line) and 
reconstructed pseudopotential (solid line) results. The reconstructed atom is at the centre and the square just encloses the 
embedding sphere. Contour levels are 0.01, 0.02, 0.025, 0.030 and 0.032 e au~ 3 , chosen to emphasise the differences between the 
original and reconstructed charge densities, (b) Valence charge density along a line in the [011] direction, with the reconstructed 
atom at the origin. Original plane-wave (dotted line) and reconstructed all-electron results (solid line) are shown. The vertical 
dashed lines are at the core radius of the pseudopotential. 



of points on the contour integral in Eq. (|22|) . The crit- 
ical source of error in the reconstruction was found to 
be the number of spherical harmonics used in the basis, 
Imax- A large enough value must be chosen for the Green 
function within the reconstruction calculation to be con- 
verged, and so to provide an accurate valence charge den- 
sity. For l max = 6 and 10, R = 0.46 % and 0.18 % re- 
spectively, and although this represents an improvement 
in the reconstruction the difference between these two re- 
sults very small. It should also be noted that the error 
introduced by a low l max is largely in the higher order 
aspherical components of the charge density - the lowest 
order terms present for aluminium (L = (00), (4m)) are 
very well converged for l max = 6. 

Second we consider the convergence of the reconstruc- 
tion with respect to the parameters of the embedding 
potential calculation. These are the number of k-points 
at which the 'band-structure' is calculated, rik; the plane- 
wave cut-off energy used in the calculations, E max ; and 
the interpolation of the spectral function applied in Eq. 
<JSJ> . In general it was found that the reconstructed charge 
density was well converged provided the embedding po- 
tential was well converged, as discussed at the beginning 
of this section. The two main points of interest are the 
convergence with respect to rik, and the sampling of the 
spectral function used to carry out the convolution inte- 
gral. 



Reconstructions were performed for the smallest 
four rik allowed within the linear analytic tetrahedron 
method, — 20, 89, 240 and 505 in the irreducible 
wedge. For these, R takes the values 2.30, 0.47, 0.44 
and 0.45 % respectively, and there is no discernable dif- 
ference between the charge densities for — 89 and 
505, suggesting that convergence has been achieved for 
89 k-points. This results is significant since the embed- 
ding potential itself was not found to be particularly well 
converged for 89 k-points. Turning to the convolution 
integral, it was found that by sampling at energy inter- 
vals of 0.3 eV below 5 eV and 2.0 eV above in Eq. ©, 
the reconstructed results are indistinguishable from those 
obtained by sampling at 0.1 eV intervals throughout the 
energy range - R takes the value 0.57 % for the more 
sparse sampling and 0.35 % for the finer sampling, and 
this difference is negligible. This coarse sampling reduces 
the number of points required from 2878 to 281, signifi- 
cantly reducing the computational effort required to con- 
struct the embedding potential. With this sampling, if, 
for example, the energy cut-off is increased from 200 eV 
to 1000 eV the computational requirements are increased 
only by a factor of ~ 3, no matter how many bands are 
included in the calculation. 

So far we have not justified using a linear interpolation 
of the spectral function in order to carry out the convolu- 
tion integral in Sec. IIIII When reconstruction results are 
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FIG. 4: Reconstructed all-electron (solid line) and original 
(dashed line, for clarity shown offset by 0.5 e Ryd~ ) den- 
sity of states within the embedding sphere. Fermi energy is 
0.207.Ryd. 

compared for a linear and cubic interpolation, and for a 
range of sampling intervals, it is found that in the use- 
ful range (ie a fine enough sampling for accurate results) 
cubic interpolation gives results only marginally better, 
or even slightly worse than linear interpolation. This 
is probably due to the spectral function not being well 
approximated by a polynomial, so a higher order inter- 
polation can give worse results than a low order method 
[30j |. It would be desirable to find a better interpolating 
function, but the diverse analytic structure due to van 
Hove singularities, and the requirement that an integral 
of the form Eq. © must be solved analytically, makes 
this a non-trivial task. 



VI. CONCLUSION 

A general method has been described for taking the 
results of a pseudopotential calculation for a given sys- 
tem and obtaining the correct all-electron charge density 
in the core region of a single atom. The reconstruction 
method provides correctly relaxed core states (they are 
not frozen to the isolated atom core states) and is in prin- 
ciple applicable to any system that can be analysed with 
plane-wave pseudopotential methods. 

The embedding potential method derived by Ingles- 
field [lT| is used, and a new analytic expression for the 
embedding potential is given here. Past applications of 
the embedding potential method have generally been lim- 
ited to models where the embedding potential is that of 



free space, or an arbitrary model embedding potential. 
This is essentially due to the difficulty of obtaining an 
accurate real space Green function for a realistic system 
- the handful of cases where a Green function (and so 
embedding potential) have been obtained from ab initio 
calculations rely on the properties of electron structure 
calculations that themselves employ a Green function (for 
example see Thisjssen and Inglesfield[4l|. Miller et al[i3|. 
Inglesfield[Hj], Cr ampin et alQ or Ishida(45|). It has 
been found that in order to construct an accurate real 
space Green function (and so embedding potential) using 
the spectral representation a complete set of eigenstates 
must be included, and here this has been taken into ac- 
count by a linear interpolation of the band structure in 
k-space and approximating the infinite number of high 
energy bands that are not available as plane-wave states. 
This results in an approximation that will converge to 
the correct form and is accurate for a realisable number 
of bands and k-points. 

Using this embedding method a localised all-electron 
atomic calculation is carried out that does not require any 
assumption of spherical symmetry in the potential neces- 
sary for previous solutions to this problem, and in essence 
makes the same physical and mathematical approxima- 
tions as state of the art all-electron density-functional 
methods such as the FLAPW method. Results of tests 
for aluminium are good - the original pseudopotential re- 
sults can be reconstructed with negligible error, and the 
all-electron results are as accurate. In addition to this 
the reconstruction accurately reproduces the density of 
states of the original substrate system. 

The major computational cost of performing a recon- 
struction lies not in the reconstruction itself, but in the 
calculation of the full set of states (ie band-structure) re- 
quired to obtain an accurate real space Green function. 
For larger systems the computational cost of obtaining 
the embedding potential from these states and perform- 
ing the reconstruction calculation are not expected to 
increase significantly, but the cost of performing the ma- 
trix diagonalisation used to obtain the full set of states 
could become prohibitive, scaling as rig. It would be 
advantageous to calculate the real space Green function 
more efficiently, and it may be possible to achieve this 
by applying iterative diagonalisation methods. Provided 
a full set of states of the self-consistent potential can be 
obtained, we expect the method presented here to be ap- 
plicable to any system. 

In a future paper [46] this method will be applied to the 
all-electron reconstruction of the core region for bulk Si, 
with the resulting charge density compared with accurate 
experimental measurements of the structure factors and 
the results of FLAPW calculations. 
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APPENDIX A: EMBEDDING METHOD 

Here we give a brief derivation of the embedding po- 
tential method of Inglesfield [ll|. This is presented to 
shed some light on the properties of the method, specif- 
ically the requirement that the embedding potential be 
an analytic function defined even where no states exist in 
the substrate system. A slightly different route is taken 
in the derivation of the embedding potential itself, which 
yields the same expressions found previously [Til \4l\ , to- 
gether with a new expression. 

1. The Embedded Hamiltonian 

Fig. [1] shows the regions I, where the new embedded 
states are obtained, and region II, the original substrate. 



Although the diagram shows the surface S as a sphere, 
this derivation also applies to any other surface. We 
proceed by finding the expectation value of the energy 
of a trial function defined in 77 as the solution of the 
Schrodinger equation for the substrate at some energy 
e, il>(r), an d in I as a trial function 4>(r). This energy 
is then expressed in a form that is dependent only on <fi 
within region I and the substrate Green function on the 
surface, S. The variational principle is then applied to 
obtain a Schrodinger equation which gives a solution in 
region / that correctly matches onto the solution in re- 
gion 77. This is a normal Schrodinger equation defined in 
/ only and with the addition of terms that are non-zero 
only at the surface S. 

The equation for ip(r) in II is 

(-V 2 + V(r) - e)V(r) = r e II, (Al) 
hence the expectation value of the energy is given by 



E = 



dip 



Jjd 3 rcj)*(j) + J n d 3 Yip*ip 
I 



(A2) 



where H is the Hamiltonian and the spatial variable is 
suppressed. The surface integral in the numerator is a 
consequence of the discontinuity in the trial function at 
S and Green's theorem, and can be interpreted as the 
contribution to the kinetic energy from this discontinuity. 
This expression has been used to apply the variational 
method to trial functions with discontinuities at a surface 
14811. but an additional condition 



(r s ) = V(r s ) 



(A3) 



is chosen here. It should be noted that no reduction in 
generality is introduced by requiring this condition to be 
satisfied since the solution in II is not explicitly described 
at any point in the derivation. 

In order to carry out a calculation localised to region 
I all explicit dependence of E on ip must be removed. 
First we remove the integral of ip over region II using 
the expression [Tlj ] 



d 3 r\ip(r) 



ii 



dh 



d dip(r s 



de dr 



(A4) 



which is an aspherical generalisation of the similar 
expression that describes the transferability of norm- 
conserving pseudopotentials [49] . This leaves the energy 
E dependent on ip only through the normal derivative at 
the surface. We express this normal derivative in terms 



of ip on the surface S as [ill ] 



dip(r s 
dn s 



[ d 2 r' s T(r s ,r' s ;e)ip(r' s ). 
Js 



(A5) 



The function r(r s ,r' s ;£) is the embedding potential, and 
by inserting Eq. (|A3IA4IA5p into Eq. (|A"2j) we obtain the 
energy of the entire system expressed entirely in terms 
of the trial function <f> in region /. The fact that an 
expression of the form Eq. (|A5|) exists is the heart of the 
embedding method. 

The condition that E is stationary with respect to 
small changes 8<p> gives the equation 



(iT + <y(r-r.) 



l l ... 



5{v-v s )j s <fv' s (T(r s ,r' s ;e) + (E-e) 



3r(r„r',;e) 
de 



E<t>{v) r el. 



(A6) 



for <f>. This takes the form of a normal Schrodinger equa- 
tion with three additional terms - the derivative term and 
the two surface integrals on the left hand side. These ex- 
tra terms take the form of a non-local potential that acts 
only at the surface S. 

In deriving Eq. (|A6[) the energy E has not been min- 
imised with respect to variations in e, hence e still ap- 
pears in the effective Schrodinger equation as a free pa- 
rameter. As a consequence of this the energy derivative 
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term is present, a first order correction to the energy at 
which the embedding potential is evaluated. This em- 
bedded Schrddinger equation will give the <fi in region / 
that is continuous with a solution in region II of energy 
e such that the combined trial function has the lowest ex- 
pectation value of energy. If the expectation value of the 
energy is further minimised with respect to s, then the 
combined trial function will be the eigenfunction of the 
complete system of eigenenergy E = e [2(|. In practice 
this requires the eigenenergy of the solution to be known 
before the equation can be directly solved, or iterative 
methods to be applied. 

Within a continuum Eq. (|A6|1 takes a simpler form as 
the energy of the state required can be chosen from the 
outset, and the first order correction is zero. In this case 
the embedded Hamiltonian for a state of energy E can 
be written 



H emb {E) = H I + S{r s -r) 



(A7) 



where H em b(E) is the embedded Hamiltonian that yields 
the states with correct boundary conditions, and Hj is 
the normal Hamiltonian for region /. The notation has 
been modified at this point such that T(r s ,r' s ; E) acting 
on a function denotes the integration over the surface iS, 
as in Eq. §M\i . 



2. Explicit Forms for the Embedding Potential 

In the above discussion the embedding potential, T, 
is defined only implicitly with no explicit form given (or 
proof that any such operator exists) . We now derive gen- 
eral forms for the embedding potential in terms of the sin- 
gle particle Green function in region //. We start with 
the defining equation for the Green function in region II, 

{-V* + V(r)-e)G(r,r';e) = 5(r-r') r, r' e II. (A8) 

Multiplying this by ?/>, and the Schrodinger equation for ip 
by G(r, r'; e), subtracting the two, integrating over region 
II and applying Green's theorem gives [llj 



i> = -G. 



dtp 
dn s 



dn' 



(A9) 



where the surface integrals are represented as matrix mul- 
tiplications. This implies a representation of the spatial 
dependence in terms of a set of basis functions which 
are orthonormal over the surface S; Q denotes this rep- 
resentation. It is straightforward to show that for this 
representation the matrix product corresponds to the in- 
tegration over the surface S. The prime indicates that 
the normal derivative of the Green function is taken with 
respect to the second spatial variable. 

The embedding potential required for the embedded 
Schrodinger equation (Eq. (|A7|) ~) is the operator that 



gives the normal derivative of a wavefunction on S in 
terms of its value on S, or 



— =T 4> 



(A10) 



where T is the matrix representation of the embedding 
potential operator. Equations (|A9|) and (|A10p lead im- 
mediately to the expression 



r = -g- 



i- 



dQ_ 

dn' 



(All) 



This general expression for T is the same as that orig- 
inally derived by Green function matching techniques 
[ill [Hol . l5lj and for von Neumann boundary conditions 
reduces to the definition given by Inglesfield in his origi- 
nal paper. 

A second expression may be obtained by taking the 
normal derivative of Eq. (|A9[) with respect to the 1 st 
spatial variable, and rearranging this to give 



r = 1 



dg_ 

dn s 



d 2 g 

dn s dn' 



(A12) 



which is an alternative and equally valid expression to Eq. 
(|A11[) for the embedding potential in terms of a Green 
function satisfying arbitrary boundary conditions on S. 
For Dirichlet boundary conditions (Q = for r or r' G 
S) this reduces to the form given by Fisher [47| . 

Finally, we note that at first glance the defining equa- 
tion (|A10[) appears to suggest that the embedding po- 
tential can be defined in terms of the eigenstates of the 
substrate system. However, this is not the case as Eq. 
(|A10p must be true for all energies in order to apply the 
variational principle, even where there are no substrate 
eigenstates. In addition to this the eigenstates of the em- 
bedded system are complex in the presence of a contin- 
uum of substrate states, hence the embedding potential 
must be an analytic function (for further discussion see 
TrailJH). 



APPENDIX B: ERROR IN THE INCOMPLETE 
SPECTRAL REPRESENTATION FOR A FREE 
ELECTRON GAS 

In order to assess the properties of the error introduced 
by lack of completeness of the spectral representation, 
the free electron case is examined. Since the high en- 
ergy Bloch states should be essentially free in character, 
it seems reasonable that this should give an indication of 
the error introduced by applying the incomplete spectral 
representation to a real system. We calculate the error as 
the contribution to the Green function from states of en- 
ergy greater than some maximum value, E cut = fcg, using 
an extended zone scheme and an expansion in spherical 
harmonics. 

The free electron Green function is spherically symmet- 
ric, this is also the case for the error, so both matrices are 
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FIG. 5: (a) Real part of free electron Green function, 5(oo,oo)i (solid line) and its normal derivative (dashed line) as a function 
of E. (b) Error in £/(oo,oo) (solid line) and its normal derivative (dashed line) at E + ie = 0.0 as a function of E cu t- 



diagonal and independent of the m index. In what fol- 
lows only one index is shown, the m index is suppressed 
and the error term is denoted £i(E, E cut ). The T matrix 
for free space is given by the imaginary part of the Green 
function for a free electron, or 



1 



Tx{E') = -kji(kr)ji(kr') 



(Bl) 



where E' = k 2 and is real [l4j. From this the error 
in the Green function, £,i(E, E cut ), can be calculated by 
applying the convolution integral (Eq. ©) to the high 
energy states to give 

jf(kr s 



&(E + ie,E cut ) = 



k 2 -(E + ie) 



dk, (B2) 



where the radial variables have been set equal to the em- 
bedding sphere radius, r s . For E cut » \E + ie\ (true for 
all E + ie we are interested in) a Taylor expansion can 
be used to give the series 



Z l (E + ie,E cut )=J2a?{E + ie) n 



where 



df '{Ecut) 



„2n-l 



kor s 



r 2n 



dx. 



(B3) 



(B4) 



These integrals can be carried out analytically by re- 
peated integration by parts, to give solutions of the form 



a\\E c 



1111 

n2n + lrf k 2 n+1 



O 



1 



(B5) 



where the second term is small for \E + ie\ <C E cut . 

Since we also require the normal derivative of the 
Green function at the surface (see Sec. [TTJ) , the error in 
this is derived in the same way. We denote the coeffi- 
cients of the normal derivative series expansion by 
and find 



-l/2r 
l/2r 2 



(B6) 



where r = r s and the last term on the right hand side is a 
step function whose value depends on which side we take 
the limit of r' — r — » in Eq. (|B1|) . This step function is 
the source of the cusp necessary for the correct analytic 
form of a Green function, and will not be present in the 
incomplete spectral representation. 

For reasonable values of E cut we find that the errors 
in both Q and J^y- are nearly constant in E + ie, so for 
convenience we examine the errors at E + ie = and 
/ = 0. Fig. [5k shows the real part of the analytic Green 
function and its normal derivative, as a function of real E 
and for I = 0. Alongside this, in Fig. [SJd, the error in the 
Green function and its normal derivative are shown as a 
function of E cut , and at E+ie = 0. These are obtained by 
evaluating a[J and 6[j analytically. It is apparent that the 
errors are significant and that convergence with respect 
to the energy cut-off is slow. For the normal derivative 
the error does not converge to zero with increasing E cut 
due to the step function in Eq. (|B6[) . For I > the errors 
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show a similar behaviour. 
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